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Abstract 



We present a new approach for Monte Carlo simulations of lattice quantum spin systems which is 
able to eliminate the negative sign problem. Its complexity is linear in the volume of the lattice. Its 
efficiency is tested on a simple 2-dimensional fermionic model. 



The investigation of quantum spin systems is important for understanding the physics 
related to strongly correlated electrons and high temperature superconductivity. A poverful 
technique to numerically study quantum spin systems is the quantum Monte Carlo (QMC) 
method (for a review see based on the Suzuki- Trotter formula Q. Unfortunately, QMC 
can suffer from the negative sign problem (NSP), which becomes exponentially serious on 
large lattices and at low temperature. In this work we present a new algorithm which is 
able to strongly reduce this problem. A more detailed description of this algorithm and a 
rigorous proof of its correctness will be published elsewhere H. This algorithm is general 
and can be applied to any quantum spin system. Its complexity is only linear in the volume 
of the lattice. We have tested its efficiency with a simple 2-dimensional fermionic model 
and we show that for this model the NSP disappears. 

We consider a quantum spin system defined on a d-dimensional finite lattice. Its direct 
simulation is not possible because the requirement in storage and work is exponential in 
the lattice size of the system. To avoid this complexity, usually one transforms it to a d+1- 
dimensional classical spin system realized by the application of the Suzuki- Trotter formula 
to the original d-dimensional system. The extra dimension is a discretization of the inverse 
temperature. We call it time. This classical spin system allows us to use a Monte Carlo 
approach for evaluating observables. The co nfigurations v defined on the d+l-dimensional 
lattice are weighted by w(v) in the partition function 



The evaluation of the weight w(v) is of small complexity, contrary to the original d- 
dimensional quantum spin system. The price to pay is that the new classical system has a 
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partition function for which generally not all weights are positive semidefinite. This funda- 
mental difficulty is usually referred to as the "negative sign" problem. It is not related to 
any approxim ations in the Monte Carlo scheme but it describes the fact that the statisti- 
cal error of the observables can become very large, increasing exponentially in the inverse 
temperature (3 and lattice volume. 



Any classical observable A is measured at the first time slice of the d+l-dimensional 
lattice by averaging over the sample of configurations generated by the Monte Carlo process. 
For simplicity, we restrict the discussion to classical observables. The ex pectation value of 
them can be written 

j: v A(v)\w(v)\-sgn{w(v)) 1 _ <Asgn> H 

<v /± > w — ^ j 7 — 7-j X / 77 — • K/'j 

}^ v \w(v)\ < sgn(w{v)) >\ w \ <sgn>\ w \ 

where < ... > w and < ... >i w i denote the average taken with the weights w and |w(i>)|, 
respectively. If the average sign is small there will be large cancelations making an accurate 
evaluation of < A > impossible so that meaningful results are dif ficult to obtain. 



We denote a Monte Carlo algorithm (j>\ w \ which generates configurations v in equilibrium 
with the distribution This algorithm can be realized in different ways. The explicit 

realization of it is for the moment irrelevant to the discussion. We denote the expectation 
value measured with the algorithm <f>t w i by eq. (Q). It is important to notice that the 
value A(v) of an observable A does not depend on the complete configuration v but only 
on the part living on the time slice at t = 0. This is a crucial property that we will use 
for constructing an algorithm whic h is able to reduce the sign problem by substantially 
increasing the average sign. 



We consider a mapping g which maps a configuration v to a configuration gv and which 
satisfies the following condition: for any observable A we have A(y) = A(gv). Such a 
mapping is easily constructed: we can require that it does not change the configuration v 
at t = 0. We consider a set Q of such mappings with the additional property that these 
mappings are bijective. It is easy to prove || that using the weight defined by 

w(v, g) = -{w(v) + w(gv)) (3) 
we obtain for the expectation value (^) of an observable 

Using this new weight for a simulation it is clear that the expectation value of A measured 
with a Monte Carlo algorithm with state space {v} and equilibrium distribution 
is the same as the one measured with a Monte Carlo algo rithm cfi^\ with state space {(v, g)} 
and equilibrium distribution \w(v,g)\. 
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We suppose that we can construct these mappings g G Q so that the average of the sign 
if measured with <^i^i satisfies 

< S3n(v) >|^| > < sgn(v) >\ w \ (5) 

where < ... >i^i and < ... >\ w \ mean the expectation taken from the \w\ and \w\ distri- 
butions, respectively. In the case that the average sign substantially increases then it is 
evident that the sign problem is suppressed. 

Such a mapping can be realized by flipping clusters of spins in v such that if w(v) ^ 
then also w(gv) 7^ and the following condition is satisfied 

condition = {(w(v) -\-w(Ecv) > 0) and 0(v,Ecv)} (6) 

where Eq denotes the flipping of the spins belonging to the cluster C. The first part of the 
condition is needed to obtain a positive weight w(v,g). The second part of the condition 
0(y, Eqv) takes the va lue true if the flipping Eq maps bijectively the configuration v to 
gv, otherwiese it takes the value false. This part of the condition is realised with very high 
accurancy using a hashing technique {|. To be spe cific, the boolean function 0(v,v') is 
defined by the following procedure 

0(v,v') : if H table (h(v')) = then H table (h(v')) := h{v) endif 
if H ta ble{h{v')) 7^ h(v) then output (D=false 
else output C=true endif 

where h is a hashing function which assigns in an arbitrary way a non vanishing integer 
label (in an interval [1, Z\ n Z) to any state v, and Htabie is a hashing table with Z integer 
entries. 

The search of the clusters follows a fixed order of operations. One can select a point in 
the d+l-dimensional lattice from a fixed list of points and construct a cluster starting from 
it. During the construction a fixed list of random numbers can be used . It is important, 
however, that the two lists always remain the same every time one applies this procedure 
to a configuration. Changing the lists is equivalent to select a new mapping in Q. If the flip 
of the constructed cluster does not satisfy th e condition (g) the next point in the list can 
be selected and a new cluster constructed and this search is repeated until the condition 
@ is satisfied. If the condition (||) is never satisfied the procedure can be stopped, when 
all the points in the list are tested once, and the original state is returned as the result. 
This algorithm defines a mapping g £ Q. The complexity of the search of a cluster is in 
this way linear in the volume of the lattice. The flip of the cluster Eq is defined so that 
A(v) = A(Ecv) for any observable A. 

This mapping defines our new weight w. A Monte Carlo algorithm 0ugi can now be used 
to update the configurations and the mappings with respect to the new weight. The first 
part of the condition @ guarantees us that after average (||) is satisfied. The second part of 
the condition (0) guarantees us that the mapping g maps a configuration v to gv bijectively 
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with probability 1 — yzj- Since the mapping is not always bijective, during the Monte 

Carlo process a systematic error is produced with probability O ■ The systematic error 
in the average (||) is then after average over the collected sample of configurations of order 
-g-. This systematic error can be tuned in the algorithm by increasing the dimensi on Z of 
the hashing table. We can chose its dimension so that is much smaller than the statistical 
error. In addition the systematic error can be measured by counting the collisions in the 
hashing table. 

If the set of mappings Q contains only one element then to avoid too much collisions 
in the hashing the dimension of the hashing table has to be much larger than the number 
of Monte Carlo iterations one desires to perform. This is doable but require s a lot of 
memory. If, however, the set of mappings Q contains a huge amount of elements then the 
hashing table can be small because the algorithm uses a selected mapping only for a short 
time and then selects a new one according to the distribution \w(v,g)\ (using a Metropolis 
acceptance test). If the set Q is big enough, the probability that the algorithm selects the 
same mapping twice is infinitesimal so that there is no need to store the history of the 
hashing tables of old mappings. 

As an example we consider a simple model of free quantum spins in a chemical potential. 
This example was analyzed by quantum Monte Carlo using the loop algorithm || in ||. We 
consider fermions living on the sites of a spatially 2-dimensional square lattice with periodic 
boundary conditions. We consider a particle with only one internal degree of freedom. 
Creation and annihilation operators c* and c x anticommute. The Hamilton operator is 

H = }(ctc x + c* , -c , ■ — etc .2 — c* , ~.c x ), (7) 

L i v x x x+i x+i x x+i x+i ■ Ln y ' 

x,i 

where i is the unit vector in the i-direction. The model is trivial and can be solved in 
momentum space. However, when one tries to simulate it with a Monte Carlo algorithm 
it shows from the algorithmic point of view many of the characteristic features of more 
complicated quantum spin systems. In addition, the exact solution of this model allows us 
to test our algorithm. We can express the grand canonical partition function 

Z = Tre~^ H -^ (8) 

(where ./V is the particle number operator) with the Suzuki- Trotter formula as explained in 

In Table 1 we present the results of the Monte Carlo simulations performed with the 
loop algorithm and with our algorithm and we compare the obtained results with the exact 
solution. The construction of the clusters in our algorithm can be performed in analogy 
with the loop algorithm ||, however using a fixed list of points in the d+l-dimensional 
lattice and a fixed list of random numbers for each mapping g € Q. For our runs we have 
used a hashing table with Z = 10000 entries. In all our simulations we have not found any 
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collisions in the hashing tables so that the systematic error on the observables is zero. 

We have applied the algorithms for various values of (5 and \x at fixed lattice spacing 
fie = 1/16. The results of both algorithms agree with the exact results within the error 
bars. It is evident that the sign problem becomes severe for the loop algorithm when the 
temperature is lowered or the chemical potential is increased. However, the sign remains 
always positive for our algorithm. 

A more exhaustive analysis of the efficiency of our algorithm is under study. Applica- 
tions of it to more physical interesting models are under way. 
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< n >o 


< n >l 


< n >e 


< sgn >o 


< sgn > l 


0.5 


2 


8 


32 


0.304(2) 


0.305(3) 


0.3049 


1.0(0) 


0.690(4) 


0.5 


4 


8 


32 


0.501(2) 


0.498(3) 


0.5000 


1.0(0) 


0.591(4) 


1.0 


2 


8 


64 


0.233(2) 


0.230(40) 


0.2321 


1.0(0) 


0.048(6) 


1.0 


3 


8 


64 


0.357(2) 


? 


0.3568 


1.0(0) 


-0.003(6) 


1.0 


4 


8 


64 


0.499(2) 


? 


0.5000 


1.0(0) 


0.002(8) 


2.0 


2 


8 


128 


0.193(2) 


? 


0.1956 


1.0(0) 


0.004(9) 


2.0 


4 


8 


128 


0.498(2) 


? 


0.5000 


1.0(0) 


-0.001(9) 


1.0 


2 


12 


64 


0.232(2) 


0.20(10) 


0.2321 


1.0(0) 


0.018(6) 


1.0 


3 


12 


64 


0.358(2) 


? 


0.3568 


1.0(0) 


0.001(3) 


1.0 


4 


12 


64 


0.501(2) 


? 


0.5000 


1.0(0) 


-0.001(4) 



Table 1: Results from the Monte Carlo simulations for different parameters. The simulations are 
done with the loop algorithm <p\ w \ and our algorithm. The expectation values of the number operator 
n and the sign obtained with the loop algorithm are denoted by < n >l and < sgn >l, the ones 
obtained with our algorithm are denoted by < n >o and < sgn >o- The mark "?" indicates that 
the measurement of the expectation value was impossible, because the statistical error exceeds the 
expectation value. The exact values of < n > are denoted by < n >e- 
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